$ontext
This code is for the following paper:
    Cai, Yongyang, William Brock, and Anastasios Xepapadeas (2022).
    Climate Change Impact on Economic Growth: Regional Climate Policy under Cooperation and Noncooperation.
    JAERE, forthcoming

Author: Yongyang Cai, The Ohio State University
Date: August 11, 2022
$offtext

parameter
    ga0N
    delaN
;

set ii /North, Tropic, South/;
set r(ii) /North, Tropic/;

* years in SSP data
set t / 2015*2515 /;

set tfirst(t);
tfirst(t)$(ord(t)=1) = YES;

$offlisting

Table Pop(t,r) 
$ondelim
$include RegionalPop_SSP1_Annual.csv
$offdelim
;

$onlisting


parameters
    gama     Capital elasticity in production function        /.300    /
    dk       Depreciation rate on capital (per year)          /.100    /

    elasmu   Elasticity of marginal utility of consumption     /1.45 /
    prstp    Initial rate of social time preference per year   /.015  /

    ga0(r)      Initial growth rate for TFP per year             
    dela(r)     Decline rate of TFP per year
    A0(r)       initial TFP
    ga(t,r)
    A(t,r)      TFP

	k0(r)
	GDP0(r)
	betas(t)   compound discount factor
;

ga0('North') = %ga0N%;
ga0('Tropic') = 0;
dela('North') = %delaN%;
dela('Tropic') = 0;

betas(t) = 1/(1+prstp)**(ord(t)-1);

* In Burke et al. (2018), GDP in 2010 is the multiplication of population in 2010 and
* the average of GDP per capita from 1980 to 2010 (after dropping those items with NA).
* However, I need to set the initial year to be 2015, so I use world bank data
* This will only affect the value of initial TFP, as we are just matching the growth rates of GDP

* Use world bank data (2015 GDP) 
GDP0('North') = 54;
GDP0('Tropic') = 19;

* RICE data in 2015
k0('North') = 100;
k0('Tropic') = 53;

ga(t,r)=ga0(r)*exp(-dela(r)*(ord(t)-1));
A0(r) = sum(tfirst, GDP0(r) / (Pop(tfirst,r)**(1-GAMA)*k0(r)**GAMA));
A(tfirst,r) = A0(r);
loop((t,r)$(ord(t)<card(t)),
  A(t+1,r) = A(t,r)/(1-ga(t,r));
);

positive variables
 K(t)            Capital stock trillions US dollars
 Y(t)
 c(t)
;
c.lo(t) = 0.01;

VARIABLES
 obj
;


EQUATIONS
 KK(t)
 GrossOutput(t)

 objfun
;

KK(t)$(ord(t)<card(t))..       K(t+1) =e= (1-DK)*K(t) + Y(t) - c(t)*Pop(t,'North');
GrossOutput(t)$(ord(t)<card(t))..	Y(t) =e= A(t,'North')*Pop(t,'North')**(1-GAMA)*K(t)**GAMA;

Objfun..            obj =e= SUM( t$(ord(t)<card(t)), betas(t)*c(t)**(1-elasmu)/(1-elasmu)*Pop(t,'North') );

K.fx(tfirst) = k0('North');

K.l(tfirst) = k0('North');
loop(t$(ord(t)<card(t)),
  Y.l(t) = A(t,'North')*Pop(t,'North')**(1-GAMA)*K.l(t)**GAMA;
  c.l(t) = 0.7*Y.l(t)/Pop(t,'North');
  K.l(t+1) = (1-DK)*K.l(t) + Y.l(t) - c.l(t)*Pop(t,'North');
);  

**********************************

** Solution options
option iterlim = 99900;
option reslim = 99999;
option solprint = off;
option limrow = 0;
option limcol = 0;
option nlp = conopt;


model Growth / all /;

solve Growth maximizing obj using nlp ;

ABORT$(Growth.MODELSTAT>2 or Growth.SOLVESTAT>1) "Model not normally completed";

parameter capY(t);
capY(t) = Y.l(t) / Pop(t,'North');

file results /output_noDam_North.csv/;
results.nd = 6 ;
results.nw = 0 ;
results.pw=20000;
results.pc=5;

put results;

loop(t$(t.val<2100),
  put capY(t):14:6; 
  put /;
);



